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ABSTRACT 

Many spacecraft attitude determination methods use exactly two vector measurements. The two vectors are typically the unit 
vector to the Sun and the Earth’s magnetic field vector for coarse “sun-mag” attitude determination or unit vectors to two 
stars tracked by two star trackers for fine attitude determination. Existing closed-form attitude estimates based on Wahba’s 
optimality criterion for two arbitrarily weighted observations are somewhat slow to evaluate. This paper presents two new 
fast quaternion attitude estimation algorithms using two vector observations, one optimal and one suboptimal. The 
suboptimal method gives the same estimate as the TRIAD algorithm, at reduced computational cost. Simulations show that 
the TRIAD estimate is almost as accurate as the optimal estimate in representative test scenarios. 

INTRODUCTION 

Suppose that we have measured two unit vectors b, and b 2 in the spacecraft body frame, e.g. unit vectors along the line of 
sight to a star or the Sun, or along the Earth’s magnetic field. We consider only unit vectors because the length of the vector 
has no information directly relevant to attitude determination. Each of these unit vectors contains two independent scalar 
pieces of attitude information. The spacecraft attitude is represented by a 3x3 proper orthogonal matrix A, i.e. A A = /, the 
3x3 identity matrix, and det A = I . Since A is thus an element of the three-parameter group SO(3), two unit vector 
measurements contain one more piece of information than is necessary to determine the attitude matrix. 

It is also necessary to know the components of the two measured vectors r j and r 2 in some reference frame. The reference 
frame is usually taken to be an inertial frame, but this is not necessary. One can use a rotating frame such as the frame 
referenced to the orbit normal vector and the local vertical. The attitude matrix to be determined is the matrix that rotates 
vectors from the reference frame to the spacecraft body frame. Thus we would like to find an attitude matrix such that 

Ar, = b, (la) 

and 

Ar 2 =b 2 . (lb) 

This is not possible in general, however, for Eq. (1) implies that 

b l b 2 =(Ar 1 )-(Ar 2 ) = r l -r 2 . (2) 

This equality is true for error-free measurements, but is not generally true in the presence of measurement errors. All 
reasonable two- vector attitude determination schemes give the same estimate when Eq. (2) is valid. 

It is clear from simple counting arguments that the two independent scalar pieces of information contained in a single vector 
measurement cannot determine the attitude uniquely. More concretely, if the attitude matrix A obeys Eq. (Ia), then so does 
the matrix fl(b P dL)Afl(r p 0 r ), for any <p b and <p r , where ff(e,0) denotes a rotation by angle <p about the axis e. This also 
makes it clear that the attitude is not uniquely determined if either the pair {b p b 2 } or the pair { r j , r 2 } is collinear. 

The earliest algorithm for determining spacecraft attitude from two vector measurements was the TRIAD algorithm 5 ' , which 
is simple to implement but does not treat the information in the two observations optimally. Many optimal estimators using 
two or more measurements have been based on a loss function introduced by Wahba 3 4 . Shuster showed a simplification ot 
his QUEST algorithm for the two-observation Wahba problem 5 , but reference 6 presented the first explicit closed-form 
solution. Optimal two-observation algorithms have not found widespread application, since they are significantly slower than 
TRIAD. In fact, algorithms specifically developed for the two-observation case may require more computations than optimal 
general purpose algorithms. Recent exceptions are Mortari’s optimal EULER-2 algorithm 7 , which approaches TRIAD in 
speed, and a suboptimal algorithm proposed by Reynolds 9 . The present paper presents two new algorithms for quaternion 
estimation from two vector measurements. The first is a very efficient optimal algorithm, which is as fast as the suboptimal 
TRIAD algorithm. The second produces the same suboptimal estimate as TRIAD, but at reduced computational cost. 



WAHB.VS PROBLE\[ 


Wahba's problem is to find the proper orthogonal matrix A that minimizes the loss function 

-Ar,|\ 


(3) 


where [b,J is a set of n unit vectors in the spacecraft body frame, { r, } are the corresponding unit vectors in the reference 
frame, and (</,} are non-negative weights. We use the invariance of the trace under cyclic permutations to rewrite Eq. (3) as 


where 


L(A) = ” trace(AZ? r ), 


(4) 

(5) 


Almost all solutions of Wahba’s problem are based on the observation that the attitude matrix that maximizes trace(AS r ) 
minimizes the loss function. The original solutions solved for the attitude matrix A directly, but most practical applications 
have been based on Davenport's ^-method 24 ' 9 , which solves for the unit attitude quaternion 10,11 
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There is a 2: l correspondence between the quaternion and the rotation matrix R(e,<p) given by 
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We will follow Shuster’s convention for quaternion products 11 , writing 
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p*q* -pq 


(9) 


This differs from the historical convention in the sign of the cross-product, and has the advantage that the order of quaternion 
multiplication is the same as the order of attitude matrix multiplication. Since the attitude matrix is a homogenous quadratic 
function of q , we can write 

\x{AB x )-q J Kq (10) 


where K is the symmetric traceless matrix 
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The ^/-method finds the optimal quaternion as the normalized eigenvector of K with the largest eigenvalue, i.e. 

Kq 0 pi “* ^jnax^opt ' 


( 12 ) 


There is no unique solution if the two largest eigenvalues of K are equal. This is not a failure of the q method; it means that 
the data aren’t sufficient to determine the attitude uniquely. 

OPTIMAL QUATERNION ESTIMATION METHOD 

In the two-observation case, it is useful to define the normalized cross products 

r. =( r i xr,)/| r | xr,| (13a) 

and 


b, h (b, x b : )/]bi x b. 


(13b) 


We note that r, or b 3 is undefined if the reference vectors or the observed vectors, respectively, are coilinear. This is the case 
noted above in which there is insufficient information to determine the attitude uniquely. It can be seen from the explicit 


solution 6 , geometrical reasoning 7 , or simply thinking about the loss function of Eq. (3) that the optimal estimate must result in 
Ar, and Ar 2 being coplanar with ^ and b 2 . This means that the optimal rotation maps the cross product vectors as 

A, P t r 3 = b 3 . (14) 

The TRIAD estimate, although not optimal, always obeys this equation; but not all two-measurement estimates do 12 . 


The quaternion rotating r 3 into b 3 with the minimum-angle rotation is, up to an overall sign, 

b 3 xr 3 ' 

1 + b 3 - r 3 _ ' 


*7 min 


J2(l + b 3 • r 3 ) 


(15) 


The most general rotation that maps r 3 into b 3 is the minimum rotation preceded by a rotation through an arbitrary angle <j> r 
about r 3 and followed by a rotation through an arbitrary angle <j> b about b 3 . As observed by Reynolds 8 , this has the quaternion 
representation 
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V2(l + b 3 -r 3 ) 


b 3 sin(0 f) /2) 
cos(0 t /2) 


b 3 x r 3 
1 + b 3 • r 3 


r 3 sin(0 r /2) 
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= <7mi„ cos(^/2) + q m sin(0/2). 
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where 0 = (fa, + 0 r and 


1 

q ' w “ V2(l + b 3 -r 3 ) 


b 3 +r 3 
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(17) 


The quaternion q ]iQ rotates the cross product vector r 3 into the body frame vector b 3 by means of a 180° rotation about the 
bisector of these two vectors. 

Now we can find the optimal quaternion by finding the angle 0 that minimizes Wahba’s loss function. Using Eqs. (4), (10), 

(16), and the half-angle formulas of trigonometry, we compute 

L(A) = a l + a 2 -(l + b 3 • r 3 ) _1 (acos0 + /?sin0), (18) 

where 

4mm^7min ” ” <?180 ^ r /l80 = (1 + ^3 ' *3) a 
and 

= c hm Kc ln m , = (1 + b 3 ■ r 3 Y l P> ( 20) 

with 

cr S (1 + b 3 • r 3 )(a,b, ■ r, + a 2 b 2 • r 2 ) + (b 3 x r 3 ) -(a,b, x r, + a 2 b 2 x r 2 ) (21) 

and 

P = (b 3 + r 3 )*(a l b 1 xr, +a 2 b 2 xr 2 ). (22) 

The loss function is minimized by setting cos 0 = a/y and sin 0 = fify , with 


(23) 

which gives 

L(A opt ) = a, + a 2 - (1 + b 3 ■ r 3 )“ 1 [ y . (24) 

The attitude quaternion estimate requires half-angles. To avoid singularities, we use cos(0/2) = -^(TTca s0) and 
sin(0/2) = sin 0/^/2(l + cos 0) when cos0>O,and sin(0/2) = -/£■(! - cos0) and cos(0/2) = sin 0/ ^2(1 - cos 0) when 
cos 0 < 0. This gives, finally, 
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for a > 0 , 


(25a) 


for a < 0. 


(25b) 



The overall sign of (he quaternion is irrelevant because of the quadratic nature of the attitude representation ot Eq. (8). This 
algorithm is very similar to Mortari’s EULER-2 solution to the two-observation Wahba problem , but it avoids explicit 
trigonometric function evaluations. 

REFERENCE FRAME ROTATIONS 


The above expressions for all the components of the optimal quaternion go to the indeterminate expression 0/0 when 
b = -r 3 . This singular condition can be avoided by solving for the attitude with respect to a reference coordinate frame 
related to the original frame by 180° rotations about the x , y, or c coordinate axis 4 * 5,13 . That is, we solve for one of the 


quaternions 


q = 


® 


-qxe, 
-q e, 


for / = L 2, 3, 


( 26 ) 


where e, is the unit vector along the I th coordinate axis. These quaternion products are trivial to implement by merely 
permuting and changing signs of the quaternion components. For example, 

q' = lq x , ft , q, , q A f ® [1, 0, 0, Of = ~ q»q 2 , - q x ] T ■ - (27) 


The equations for the inverse transformations are the same, since a 180° rotation in the opposite direction has the same effect. 
These rotations are also easy to implement on the input data, since a rotation about axis i simply changes the signs of the 7 th 
and k m columns of r,, r 2 , and r 3 , where { i, y, &} is a permutation of { 1, 2, 3}. 

The original QUEST implementation performed sequential rotations one axis at a time, until an acceptable reference 
coordinate system was found. It is clearly preferable to save computations by choosing a single desirable rotation as early in 
the computation as possible. Consider the effect on the inner product b 3 r 3 of rotating the reference frame about the ** axis 


(b 3 ‘ Ts ) routed (^3 (*3 Junrotated [ 2 (b 3 ), ( **3 ) / ^ 3 r 3 3 un rouicU 

We can maximize this inner product in the rotated frame by performing no rotation if b 3 * r 3 is the maximum of 
{ (b 3 ) r (r 3 )., (b 3 ) y (r 3 ) j| ,(b 3 ) t (r 3 ) t , b 3 r 3 } xuuolatedi while a rotation about the 1 th axis is performed if (b 3 ) i (r 3 ) / is the maximum. 
This will ensure the largest value for b 3 • r 3 in the rotated frame. The rotation is easily “undone” by Eq. (27) or its equivalent 
after the quaternion has been computed. 

SUBOPTIMAL QUATERNION ESTIMATION METHOD 

In the limit that the first measurement is far more accurate than the second, Eq. (la) is satisfied exactly. This is often the case 
of interest, and is the case generally treated by the TRIAD algorithm 1 ' 2 . The quaternion estimate taking r, into b} is given, in 
analogy with Eq. (16), by 

b| + r t 


(28) 
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cos(v^/2) 


b, xr, 

1 + b, r 


+ sin(v/"/2) 
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(29) 


We compute the TRIAD-equivalent quaternion by finding the angle t// - that minimizes Wahba’s loss function. In parallel with 
the optimal case, we find 

L{A) = a,[l -(b, ■ b 2 )(r, • r, ) — (I -t- b, • r^ipcosy/ + vsin y/)], (30) 

where 

H = (1 + b, r,)[(b, xbj)-(r, xr,)]-[b, (r, xr,)][r, (b, xb,)] (31) 

and 

v = (b, +r,) [(b, xb 2 )x(r, xr,)]. (32) 

The loss function is minimized by setting cost// = jj.jp and sint// = v/p. with 


p = 7 m 2 + v 2 , 


(33) 


which gives 


a.V, AD ) = a.[l -(b, b,)(r, r,)-(l + b, r, )“'p ]. 


(34) 


The attitude quaternion estimate is 


and 


^/triad “ 


1 pp + jUHb, xr,) + V(b, - r, ) 

2>(p + ^)(l + b, r,)! (P + P)(l + b, r,) 


for u > 0 , 


(35a) 


*7triad “ 


2^/p(p-p)(l + b, r,) 


v(b, x r, ) + (p - p)(b, + r, ) 
v(I + b,-r,) 


for ju < 0 . 


(35b) 


This estimate is less expensive to compute than the optimal estimate, since it is not necessary to normalize the cross products 
r t x r, and bj x b 2 . However, Eqs. (31) and (32) show that both jx and vare zero and the estimate is undefined if either of 
these cross products is zero, indicating collinearity of the reference or body frame vectors. The case of =-r, is handled by 
reference frame rotations in the same manner as for the optimal estimate. 


TESTS 


MATLAB implementations of these algorithms were coded and tested. The computational speed of the new optimal 
algorithm is 156 or 158 floating point operations, depending on whether or not a reference frame rotation is required. These 
are exactly the same number of flops as used by the non-optimal TRIAD algorithm followed by the extraction of a quaternion 
from the TRIAD attitude matrix. The effort for the optimal algorithm includes the computation of the loss function, which is 
not provided by TRIAD, however. The suboptimal quaternion estimation algorithm required 107 or 109 flops, without 
computation of the loss function. 

For accuracy tests, we simulate 1000 test cases with uniformly distributed random attitude quaternions. For each attitude we 
generate two random observation vectors, independently and uniformly distributed on the unit sphere. We use inverse of the 
attitude quaternion to map the observation vectors to the reference frame. The reference vectors are corrupted by Gaussian 
random noise with specified standard deviations and then normalized. The optimal estimation algorithm uses the inverse 
squared measurement standard deviations as the weights in Wahba’s loss function. 

Two different test scenarios were simulated. The first scenario has measurement standard deviations of 1 arc minute per axis 
on the first vector and 2° on the second, which would be appropriate for measurements of a fine digital sun sensor and a 
triaxial magnetometer, respectively. The attitude estimation errors for this scenario, the rotation angles between true attitude 
and the estimates, are plotted in Figure 1 as a function of [b, xb 2 |. Since the errors are expected to be roughly inversely 
proportional to \b { x b 2 |, the errors have been multiplied by this quantity for plotting, allowing the errors to be plotted on a 
more-or-less uniform scale for the entire range of angles between the observed vectors 6 . The points for the suboptimal 
estimation errors and the optimal errors fall on top of each other in Figure 1; they are indistinguishable at the plotting 
accuracy. This was to be expected with the large the ratio of the measurement weights, 120 2 , in this scenario. 

The second scenario, shown in Figure 2, used equal measurement standard deviations of 2° on both, as would be appropriate 
for measurements of a coarse sun sensor and a triaxial magnetometer. This scenario clearly exhibits differences between the 
suboptimal optimal estimation errors. This is not surprising, since the difference between the suboptimal TRIAD estimator 
and the optimal is greatest in the case of equal weights. Interestingly enough, the errors agree very closely when they are 
greater than 20°, as shown by the points above the diagonal line from (0, 0) to (0.7, 14) in Figure 2. These solutions with 
large errors are not very useful, however. In fact, the figures show that two-observation estimators should not be used for 
|b t x b,| < 0.4 if it is desired to limit the estimation errors to be less than 10° in the first scenario or in 20° the second. 

It appears that the optimal estimates have smaller errors than the TRIAD estimates in the second scenario, but this is difficult 
to quantify by examining Figure 2. To establish the advantage of the optimal estimator, the probability distribution of 
|bj x b,| times the estimation errors in 10,000 simulations of the second scenario are plotted in Figure 3. It is clear from this 
figure that the optimal estimator offers only marginal gains over the TRIAD estimate. At the 95% confidence level, |b L x b 2 | 
times the errors of the optimal and suboptimal algorithms are less than 5.3° and 5.6°, respectively; these scaled errors are less 
than 6.7° and 6.9° at the 99% confidence level. 



Figure 1. Estimation errors in scenario with unequal measurement errors 


CONCLUSIONS 


We have found two new, fast quaternion estimation methods using exactly two vector measurements. These methods are 
applicable to a variety of problems, including coarse “sun-mag” attitude estimation using the unit vector to the Sun and the 
Earth’s magnetic field vector and precise estimation using unit vectors to stars tracked by two star trackers. The first of the 
algorithms is optimal in Wahba’s sense and is as fast as the TRIAD algorithm, even including computation of Wahba’s loss 
function, which TRIAD does not provide. The second, suboptimal algorithm provides estimates identical to TRIAD with 
reduced computational effort. The accuracy of these algorithms w r as examined in simulations of two scenarios, one with 
equal measurement errors and one with errors differing by two orders of magnitude. These simulations show that the TRIAD 
estimate is almost as accurate as the optimal estimate in both scenarios. 



Figure 2. Estimation errors in scenario with equal measurement errors 
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Left curve for optimal estimator, right curve for suboptimal. 
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